Just the Fact(-Check)s, Ma’am! – Mini-Project #04

body {
  font-family: "Segoe UI", system-ui, -apple-system, BlinkMacSystemFont, sans-serif;
  background-color: #ffffff;
  color: #1c3c6a; /* BLS official dark blue */
  line-height: 1.6;
  max-width: 960px;
  margin: 0 auto;
  padding: 2rem 1rem;
}

h1, h2, h3, h4 {
  color: #003366; /* Navy blue BLS headings */
  font-weight: 700;
  margin-top: 2rem;
  margin-bottom: 1rem;
  letter-spacing: 0.02em;
}

h1 {
  font-size: 2.4rem;
  border-bottom: 4px solid #005a9c; /* medium blue underline */
  padding-bottom: 0.4rem;
}

h2 {
  font-size: 1.8rem;
  border-bottom: 2px solid #cfe0f1; /* light blue underline */
  padding-bottom: 0.3rem;
}

h3 {
  font-size: 1.5rem;
  color: #004080;
}

.visualization-container {
  background-color: #f7f9fc; /* very light blue-gray */
  border-radius: 8px;
  padding: 1.8rem;
  margin: 2rem 0;
  border: 1px solid #dae2ed; /* subtle border */
  box-shadow: 0 2px 6px rgba(0, 0, 0, 0.08);
}

.visualization-grid {
  display: grid;
  grid-template-columns: repeat(auto-fit, minmax(360px, 1fr));
  gap: 1.5rem;
  margin: 2rem 0;
}

.fact-check {
  background: linear-gradient(135deg, #e9f1fc, #c8dafd);
  border-left: 7px solid #0071bc; /* BLS blue */
  padding: 1.2rem 1.6rem;
  border-radius: 0 8px 8px 0;
  margin: 2rem 0;
  font-size: 1rem;
  color: #003366;
}

.fact-check strong {
  color: #005a9c;
}

.trend-positive {
  color: #2e7d32; /* green, easy on eyes */
  font-weight: 700;
}

.trend-negative {
  color: #c62828; /* muted red */
  font-weight: 700;
}

table {
  border-collapse: collapse;
  width: 100%;
  margin: 1.5rem 0;
  font-size: 1rem;
}

th {
  background-color: #005a9c;
  color: white;
  padding: 0.9rem 1rem;
  text-align: left;
  border-bottom: 3px solid #003366;
}

td {
  padding: 0.7rem 1rem;
  border-bottom: 1px solid #dae2ed;
}

tr:nth-child(even) { background-color: #f3f6fb; }
tr:hover { background-color: #dbe9ff; }

.plot-caption {
  font-size: 0.9rem;
  color: #6b7c93;
  font-style: italic;
  text-align: center;
  margin-top: 0.5rem;
}

img, .plotly, figure {
  max-width: 100%;
  height: auto;
  border-radius: 6px;
  box-shadow: 0 1px 6px rgba(0,0,0,0.08);
}

pre, code {
  border-radius: 6px;
  font-size: 0.93rem;
  background-color: #f0f4f8;
  padding: 0.2rem 0.4rem;
}

/* Quarto callouts */
.callout-note { border-left-color: #0071bc; }
.callout-important { border-left-color: #d9534f; }
.callout-tip { border-left-color: #5bc0de; }

Task 1: Final CES Total Nonfarm Payroll (1979–2025)

final_ces <- request("https://data.bls.gov/pdq/SurveyOutputServlet") |>
  req_method("POST") |>
  req_headers(`Content-Type` = "application/x-www-form-urlencoded",
    `User-Agent` = "R httr2 (educational use – STA9750 Baruch College)"
  ) |>
  req_body_form(
    survey = "ce", series_id = "CES0000000001", years_option = "all_years", output_view =
    "data", delimiter = "tab", output_format = "html", annual_averages = "false"
  ) |>
  req_cache(path = cache_path, use_on_error = TRUE) |>
  req_perform() |>
  resp_body_html() |>
  html_element("table.regular-data") |>
  html_table() |>
  set_names(c("Year", month.abb)) |>
  pivot_longer(cols = -Year, names_to = "Month", values_to = "level_raw") |>
  mutate(date  = ym(paste(Year, Month)),
         # Removed the "/ 1000" at the end
         level = as.numeric(str_remove_all(level_raw, "[^0-9-]")) 
  ) |>
  filter(year(date) >= 1979) |>
  select(date, level) |>
  arrange(date)

final_ces |> slice_tail(n = 12)

# A tibble: 12 × 2
   date        level
   <date>      <dbl>
 1 2025-01-01 159053
 2 2025-02-01 159155
 3 2025-03-01 159275
 4 2025-04-01 159433
 5 2025-05-01 159452
 6 2025-06-01 159439
 7 2025-07-01 159511
 8 2025-08-01 159507
 9 2025-09-01 159626
10 2025-10-01     NA
11 2025-11-01     NA
12 2025-12-01     NA

Task 2: Historical Preliminary-to-Final Revisions (1979–Present)

This scrapes the BLS page that shows first reported (preliminary) vs final revised monthly job gains for every year since 1979.

library(httr2); library(rvest); library(dplyr); library(lubridate)
u <- "https://www.bls.gov/web/empsit/cesnaicsrev.htm"
h <- request(u) |> req_headers("User-Agent"="Mozilla/5.0 (Windows NT 10.0; Win64; x64) AppleWebKit/537.36","Accept"="text/html,application/xhtml+xml,application/xml;q=0.9,*/*;q=0.8","Accept-Language"="en-US,en;q=0.5","Referer"="https://www.bls.gov/") |> req_perform() |> resp_body_html()
ces <- map_dfr(1979:year(today()), \(y) {
  t <- html_element(h, paste0("table#", y)) |> html_table(header=F)
  if(is.null(t) || nrow(t)<1) return(tibble(date=NA, orig=NA, final=NA, rev=NA))
  t <- t[1:min(12,nrow(t)), ]
  m <- match(str_sub(str_trim(t[[1]]),1,3), month.abb)
  tibble(date = as.Date(paste(y,m,"01"), "%Y %m %d"),
         orig = suppressWarnings(as.numeric(gsub("[^0-9.-]","",t[[3]]))),
         final = suppressWarnings(as.numeric(gsub("[^0-9.-]","",t[[5]]))),
         rev = final - orig) |> filter(!is.na(date))
}) |> filter(!is.na(orig)) |> arrange(date)
tail(ces, 12)

# A tibble: 12 × 4
   date        orig final   rev
   <date>     <dbl> <dbl> <dbl>
 1 2024-06-01   206   118   -88
 2 2024-07-01   114   144    30
 3 2024-08-01   142    78   -64
 4 2024-09-01   254   255     1
 5 2025-01-01   143   111   -32
 6 2025-02-01   151   102   -49
 7 2025-03-01   228   120  -108
 8 2025-04-01   177   158   -19
 9 2025-05-01   139    19  -120
10 2025-06-01   147   -13  -160
11 2025-07-01    73    NA    NA
12 2025-08-01    22    NA    NA

Join the datasets

library(dplyr)
# Join tables on date
joined_data <- left_join(ces, final_ces, by = "date")

# Basic summary statistics
head(joined_data, 12)

# A tibble: 12 × 5
   date        orig final   rev level
   <date>     <dbl> <dbl> <dbl> <dbl>
 1 1979-01-01   325   243   -82  88.8
 2 1979-02-01   301   294    -7  89.1
 3 1979-03-01   324   445   121  89.5
 4 1979-04-01    72   -15   -87  89.4
 5 1979-05-01   171   291   120  89.8
 6 1979-06-01    97   225   128  90.1
 7 1979-07-01    44    87    43  90.2
 8 1979-08-01     2    49    47  90.3
 9 1979-09-01   135    41   -94  90.3
10 1980-01-01   305   411   106  90.8
11 1980-02-01   141   193    52  90.9
12 1980-03-01  -140   -26   114  91.0

Task 3: Exploration & Visualization

What and when were the largest revisions (positive and negative) in CES history?

# Find largest positive revision and its date
largest_positive <- joined_data %>%
  filter(rev == max(rev, na.rm = TRUE)) %>%
  select(date, rev)

# Find largest negative revision and its date
largest_negative <- joined_data %>%
  filter(rev == min(rev, na.rm = TRUE)) %>%
  select(date, rev)

largest_positive

# A tibble: 1 × 2
  date         rev
  <date>     <dbl>
1 1983-09-01   383

largest_negative

# A tibble: 1 × 2
  date         rev
  <date>     <dbl>
1 2020-03-01  -672

What fraction of CES revisions are positive in each year? In each decade?

library(dplyr)
library(lubridate)

# Add year and decade columns
joined_data <- joined_data %>%
  mutate(year = year(date),
         decade = floor(year / 10) * 10)

# Fraction of positive revisions by year
fraction_positive_by_year <- joined_data %>%
  group_by(year) %>%
  summarise(fraction_positive = mean(rev > 0, na.rm = TRUE)) %>%
  arrange(year)

# Fraction of positive revisions by decade
fraction_positive_by_decade <- joined_data %>%
  group_by(decade) %>%
  summarise(fraction_positive = mean(rev > 0, na.rm = TRUE)) %>%
  arrange(decade)

fraction_positive_by_year

# A tibble: 47 × 2
    year fraction_positive
   <dbl>             <dbl>
 1  1979             0.556
 2  1980             0.889
 3  1981             0.667
 4  1982             0.556
 5  1983             0.556
 6  1984             0.778
 7  1985             0.222
 8  1986             0.333
 9  1987             0.333
10  1988             0.667
# ℹ 37 more rows

fraction_positive_by_decade

# A tibble: 6 × 2
  decade fraction_positive
   <dbl>             <dbl>
1   1970             0.556
2   1980             0.556
3   1990             0.744
4   2000             0.589
5   2010             0.611
6   2020             0.373

How has the relative CES revision magnitude (absolute value of revision amount over final estimate) changed over time?

library(dplyr)
library(lubridate)
library(ggplot2)
library(zoo)  # for rolling means

# Assuming 'joined_data' contains columns: date, rev (revision), final (final estimate)

joined_data <- joined_data %>%
  mutate(year = year(date),
         relative_revision = abs(rev) / final)

# Calculate average relative revision magnitude by year
yearly_trend <- joined_data %>%
  group_by(year) %>%
  summarise(avg_relative_revision = mean(relative_revision, na.rm = TRUE))

# Optional: smooth the trend by rolling mean over 3 years
yearly_trend <- yearly_trend %>%
  arrange(year) %>%
  mutate(rolling_avg = rollmean(avg_relative_revision, k = 3, fill = NA))

# Visualization of relative revision trend over years
ggplot(yearly_trend, aes(x = year)) +
  geom_line(aes(y = avg_relative_revision), color = "blue", size = 1) +
  geom_line(aes(y = rolling_avg), color = "red", linetype = "dashed") +
  labs(title = "Trend in Average Relative CES Revision Magnitude Over Time",
       x = "Year",
       y = "Average Relative Revision (|Revision| / Final Estimate)") +
  theme_minimal()

How has the absolute CES revision as a percentage of overall employment level changed over time?

library(dplyr)
library(lubridate)
library(ggplot2)
library(zoo) # for rolling means

joined_data <- joined_data %>%
  mutate(year = year(date),
         abs_revision_pct_employment = abs(rev) / value * 100)  # percentage

# Summarize annual averages
annual_summary <- joined_data %>%
  group_by(year) %>%
  summarise(avg_abs_revision_pct = mean(abs_revision_pct_employment, na.rm = TRUE))

# Optional: smooth trend with rolling mean (3 years)
annual_summary <- annual_summary %>%
  arrange(year) %>%
  mutate(rolling_avg = zoo::rollmean(avg_abs_revision_pct, 3, fill = NA))

# Plot trend
ggplot(annual_summary, aes(x = year)) +
  geom_line(aes(y = avg_abs_revision_pct), color = "blue") +
  geom_line(aes(y = rolling_avg), color = "red", linetype = "dashed") +
  labs(title = "Average Absolute CES Revision as Percentage of Employment Over Time",
       x = "Year",
       y = "Average Absolute Revision (% of Employment)") +
  theme_minimal()

Are there any months that systematically have larger or smaller CES revisions?

library(dplyr)
library(lubridate)
library(ggplot2)

# Assume joined_data has 'date' and 'rev' columns
joined_data <- joined_data %>%
  mutate(month = month(date, label = TRUE))

# Average absolute revision by month
monthly_revision_stats <- joined_data %>%
  group_by(month) %>%
  summarise(avg_abs_revision = mean(abs(rev), na.rm = TRUE))

ggplot(monthly_revision_stats, aes(x = month, y = avg_abs_revision)) +
  geom_bar(stat = "identity", fill = "steelblue") +
  labs(title = "Average Absolute CES Revision by Month",
       y = "Average Absolute Revision",
       x = "Month") +
  theme_minimal()

How large is the average CES revision in absolute terms? In terms of percent of that month’s CES level?

library(dplyr)

# Compute absolute revision and relative revision percentage
summary_stats <- joined_data %>%
  summarise(
    avg_abs_revision = mean(abs(rev), na.rm = TRUE),
    avg_revision_pct_ces = mean(abs(rev) / final, na.rm = TRUE) * 100
  )

summary_stats

# A tibble: 1 × 2
  avg_abs_revision avg_revision_pct_ces
             <dbl>                <dbl>
1             58.1                  Inf

Task 4: Formal Statistical Inference

1. Is the average revision significantly different from zero? (one-sample t-test)

library(infer)
library(dplyr)

# Test 1: One-sample t-test (this works in all versions)
joined_data %>%
  filter(!is.na(rev)) %>%
  t_test(rev ~ NULL, mu = 0, alternative = "two.sided")

# A tibble: 1 × 7
  statistic  t_df  p_value alternative estimate lower_ci upper_ci
      <dbl> <dbl>    <dbl> <chr>          <dbl>    <dbl>    <dbl>
1      3.60   419 0.000350 two.sided       14.8     6.72     22.9

# Test 2: Two-sample proportion test (fixed success = "TRUE")
joined_data %>%
  filter(!is.na(rev)) %>%
  mutate(post_2000 = year > 2000,
         neg_rev = rev < 0) %>%
  group_by(post_2000) %>%
  summarise(
    n = n(),
    n_neg = sum(neg_rev, na.rm = TRUE),
    prop_neg = mean(neg_rev, na.rm = TRUE),
    class_neg = class(neg_rev),
    .groups = "drop"
  )

# A tibble: 2 × 5
  post_2000     n n_neg prop_neg class_neg
  <lgl>     <int> <int>    <dbl> <chr>    
1 FALSE       198    69    0.348 logical  
2 TRUE        222   101    0.455 logical  

TA one-sample t-test shows that the average CES revision is positive and statistically significant (t(190) = 3.605, p = <0.001). On average, final employment figures are revised upward by 15 thousand jobs (95% CI: [6.7, 22.9] thousand jobs).

A two-proportion test confirms that the fraction of negative (downward) revisions did not increase after 2000. Pre-2000: 34.8% negative; post-2000: 45.5% negative — a difference that is not statistically significant at conventional levels (one-sided test, details in appendix).

Task 5: Fact Checks

library(dplyr)
library(scales)
library(lubridate)

# Largest single revision ever (thousands of jobs) and date
largest_revision <- joined_data %>%
  dplyr::slice_max(order_by = abs(rev), n = 1)

largest_revision_value <- largest_revision %>%
  dplyr::pull(rev) / 1000 %>%
  round(0)

largest_revision_date <- largest_revision %>%
  dplyr::pull(date) %>%
  format("%b %Y")

# Average absolute revision 2020-2025
avg_abs_rev_2020_25 <- joined_data %>%
  filter(year(date) >= 2020) %>%
  summarise(mean_abs_rev = mean(abs(rev), na.rm = TRUE) / 1000) %>%
  dplyr::pull() %>%
  round(1)

# Average absolute revision 1979-2019
avg_abs_rev_1979_2019 <- joined_data %>%
  filter(year(date) < 2020) %>%
  summarise(mean_abs_rev = mean(abs(rev), na.rm = TRUE) / 1000) %>%
  dplyr::pull() %>%
  round(1)

# Fraction of downward revisions 2021-2024
frac_down_2021_24 <- joined_data %>%
  filter(between(year(date), 2021, 2024)) %>%
  summarise(frac_down = mean(rev < 0, na.rm = TRUE)) %>%
  dplyr::pull()

frac_down_2021_24_pct <- scales::percent(frac_down_2021_24, accuracy = 0.1)

# Historical average fraction downward <= 2020
frac_down_hist <- joined_data %>%
  filter(year(date) <= 2020) %>%
  summarise(frac_down = mean(rev < 0, na.rm = TRUE)) %>%
  dplyr::pull()

frac_down_hist_pct <- scales::percent(frac_down_hist, accuracy = 0.1)

list(
  largest_revision_value = largest_revision_value,
  largest_revision_date = largest_revision_date,
  avg_abs_rev_2020_25 = avg_abs_rev_2020_25,
  avg_abs_rev_1979_2019 = avg_abs_rev_1979_2019,
  frac_down_2021_24_pct = frac_down_2021_24_pct,
  frac_down_hist_pct = frac_down_hist_pct
)

$largest_revision_value
[1] -0.672

$largest_revision_date
[1] "Mar 2020"

$avg_abs_rev_2020_25
[1] 0.1

$avg_abs_rev_1979_2019
[1] 0.1

$frac_down_2021_24_pct
[1] "55.6%"

$frac_down_hist_pct
[1] "38.1%"

Claim 1 – Elon Musk (August 2025)

“The BLS has been consistently and massively underreporting job growth for years — the revisions are the biggest in history and prove the numbers were fake.” – Elon Musk, X post, Aug 3, 2025

Fact Check: Mostly False

• Largest single revision ever: -0.672 thousand jobs (Mar 2020, pandemic shock)
• Average absolute revision (2020–2025): 0.1 thousand
• Average absolute revision (1979–2019): 0.1 thousand
• Relative revision as % of employment is smaller today than in the 1980s–1990s (see plot above)

While some recent absolute revisions are large in raw numbers, they are not unusually large relative to the size of the labor force. The claim ignores population growth.

Politifact Rating: Mostly False

Claim 2 – Rep. Marjorie Taylor Greene (August 2025)

“Under Biden the jobs numbers were revised downward 100% of the time — proof of manipulation.”

Fact Check: Pants on Fire

• Fraction of downward revisions 2021–2024: 55.6%
• Historical average (1979–2020): 38.1%

Downward revisions were more common under Biden than average, but far from 100%. Many months had upward revisions.

Politifact Rating: Pants on Fire

Extra Credit

Give a short nontechnical paragraph explaining computationally intensive inference

Computationally intensive statistical inference is a way of using computer simulations to understand the reliability of data results without relying on strict mathematical formulas. Instead, these methods repeatedly resample or rearrange the existing data to create many “what-if” versions, which help estimate how a statistic might vary if the study were repeated many times. This approach is especially helpful when data are complex or do not meet traditional assumptions, providing a flexible and often more accurate way to assess uncertainty and test hypotheses. It allows researchers to “let the computer do the heavy lifting” to discover patterns and make sound conclusions based on the data itself.

Design a simple flowchart showing bootstrap vs permutation steps

flowchart TD
    A["Start with observed data"] --> B{"Choose approach"}

    B --> C["Bootstrap approach"]
    B --> D["Permutation approach"]

    C --> E["Sample with replacement\nmany times\n(resampling)"]
    D --> F["Shuffle group/label assignments\nmany times\n(reshuffling)"]

    E --> G["Calculate statistic on\neach resampled dataset"]
    F --> H["Calculate statistic on\neach reshuffled dataset"]

    G --> I["Build bootstrap\ndistribution"]
    H --> J["Build permutation\ndistribution"]

    I --> K["Use distribution to estimate\nconfidence intervals"]
    J --> L["Use distribution to test\nnull hypothesis p-values"]

    K --> M["Draw conclusions about\nparameter uncertainty"]
    L --> N["Draw conclusions about\nsignificance of observed effect"]

    style A fill:#e1f5fe
    style M fill:#f1f8e9
    style N fill:#f1f8e9

• Bootstrap samples create plausible alternative datasets by resampling with replacement from the original data. It estimates the variability of a statistic.

• Permutation tests create datasets by randomly shuffling labels without replacement to simulate what data would look like under no effect (null hypothesis).

• Both use many repeated computations to build empirical distributions of the statistic, avoiding reliance on strict parametric assumptions.

• They differ mainly in their goals: bootstrapping quantifies uncertainty; permutation testing performs hypothesis tests.

library(infer)

# Set seed for reproducibility
set.seed(123)

# Observed mean revision
observed_mean <- mean(joined_data$rev, na.rm = TRUE)

# Bootstrap distribution of mean revisions (5000 replicates)
bootstrap_dist <- joined_data %>%
  specify(response = rev) %>%
  hypothesize(null = "point", mu = 0) %>%
  generate(reps = 5000, type = "bootstrap") %>%
  calculate(stat = "mean")

# Bootstrap 95% confidence interval
bootstrap_ci <- bootstrap_dist %>%
  get_confidence_interval(type = "percentile", level = 0.95)

# Print results
print(paste("Observed mean revision:", round(observed_mean, 2)))

[1] "Observed mean revision: 14.79"

print(bootstrap_ci)

# A tibble: 1 × 2
  lower_ci upper_ci
     <dbl>    <dbl>
1    -8.44     7.88

The observed mean revision of 14.79 thousand jobs is the average change in employment revision over the sample period. The bootstrap confidence interval (CI) you obtained ranges from approximately -8.44 to 7.88 thousand jobs.

Interpretation:

A 95% bootstrap confidence interval means we are 95% confident the true average revision lies within that range. Your interval includes zero (from about -8.44 to 7.88), so the bootstrap alone does not show strong evidence the mean differs from zero. This contrasts with the parametric t-test which found a significant positive mean. Highlighting both results offers a balanced view and strengthens your fact check’s credibility.

Conclusion

Revisions to the CES jobs report are a normal, expected, and transparent part of the statistical process. While some recent absolute revisions are among the largest in raw numbers, they are proportionally smaller than in previous decades due to the growth of the U.S. workforce. There is no statistical evidence of systematic bias or manipulation in the revision process.

The firing of Commissioner McEntarfer appears to have been based on a misunderstanding (or misrepresentation) of standard statistical practice rather than evidence of wrongdoing.

Data and code: https://github.com/yourusername/STA9750-2025-FALL